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A many-body theory of molecular junction transport based on nonequilibrium Green's functions 
is developed, which treats coherent quantum effects and Coulomb interactions on an equal footing. 
The central quantity of the many-body theory is the Coulomb self-energy matrix Sc of the junction. 
Ec is evaluated exactly in the sequential tunneling limit, and the correction due to finite tunnel- 
ing width is evaluated self-consistently using a conserving approximation based on diagrammatic 
perturbation theory on the Keldysh contour. Our approach reproduces the key features of both 
the Coulomb blockade and coherent transport regimes simultaneously in a single unified transport 
theory. As a first application of our theory, we have calculated the thermoelectric power and dif- 
ferential conductance spectrum of a benzenedithiol-gold junction using a semi-empirical 7r-electron 
Hamiltonian that accurately describes the full spectrum of electronic excitations of the molecule up 
to 8-10eV. 

PACS numbers: 73.63.-b, 85.65.+h, 72.10.Bg, 31.15.Nc 



I. INTRODUCTION 

Electron transport in single-molecule junctions 1 -! 2 -^ is of fundamental interest as a paradigm for nanosystems far 
from equilibrium, and as a means to probe important chemical^ and biological^ processes, with myriad potential device 
applications^^ A general theoretical framework to treat the many-body problem of a single molecule coupled to 
metallic electrodes does not currently exist. Mean-field approaches based on density-functional theor y 8 ! 9 : 10 ! 11 : 12 ! 13 : 14 — 
the dominant paradigm in quantum chemistry — have serious shortcoming s 15 ! 16 ' 17 because they do not account for 
important interaction effects like Coulomb blockade. 

An alternative approach is to solve the few-body molecular Hamiltonian exactly, and treat electron hopping between 
molecule and electrodes as a perturbation. This approach has been used to describe molecular junction transport in 
the sequential-tunneling regime ; 15 ' 18 ! 19 but describing coherent quantum transport in this framework remains an open 
theoretical problem. Higher-order tunneling processes may be treated rigorously in the density-matrix formalism ; 20 ' 21 
but the expansion is typically truncated at fourth order ; 20 ' 22 the calculation of higher-order terms being prohibitively 
difficult. 

In this article, we develop a many-body theory of molecular junction transport based on nonequilibrium Green's 
function s 23 ' 24 ' 25 (NEGF), in order to utilize physically motivated approximations that sum terms of all orders. The 
junction Green's functions are calculated exactly in the sequential-tunneling limit, and the corrections to the electron 
self-energy due to finite tunneling width are included via Dyson-Keldysh equations. The tunneling self-energy is 
calculated exactly using the equations-of-motions method ; 2 ' 26 while the correction to the Coulomb self-energy is 
calculated using diagrammatic perturbation theory. In this way, tunneling processes are included to infinite order, 
meaning that any approximation utilized is a truncation in the physical processes considered rather than in the order 
of those processes. 

Our approach reproduces the key features of both the Coulomb blockade and coherent transport regimes simultane- 
ously in a single unified transport theory. Nonperturbative effects of intramolecular correlations are included, which 
are necessary to accurately describe the HOMO-LUMO gap, essential for a quantitative theory of transport. 

As a first application of our many-body transport theory, we investigate the benchmark system of benzene(l,4)dithiol 
(BDT) with gold electrodes. Two key parameters determining the lead-molecule coupling — the tunneling width T 
and the chemical potential offset A/x — are fixed by comparison to linear-response measurements of the thermoelec- 
tric powe r 27 ' 28 and electrical conductance^ The nonlinear junction response is then calculated. The differential 
conductance spectrum of the junction exhibits an irregular "molecular diamond" structure analogous to the regular 
Coulomb diamonds 30 * observed in quantum dot transport experiments, as well as clear signatures of coherent quantum 
transport — such as transmission nodes due to destructive interference — and of resonant tunneling through molecular 
excited states. 

The article is organized as follows: A detailed derivation of the many-body theory of molecular junction transport is 
presented in Sec. [TXJ Two useful approximate solutions for the nonequilibrium Coulomb self-energy are given in Sees. 
HID II and HID 21 The details of the model used to describe a molecular heterojunction consisting of a 7r-conjugated 
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molecule covalently bonded to metallic electrodes are presented in Sec. IIIII The electric and thermoelectric response 
of a BDT-Au junction are calculated in Sec. HVi and compared to experimental results. A discussion and conclusions 
are presented in Sec. [Vj 



II. NONEQUILIBRIUM MANY-BODY FORMALISM 
A. Molecular junction Hamiltonian 

The Hamiltonian of a junction consisting of a molecule coupled to several metallic electrodes may be written 

Hj unction = i^mol + H leads + Ht- (1) 

The molecular Hamiltonian can be formally divided into one-body and two-body terms -ff mo i = H^Jl + H^h. In 
general, neglecting spin-orbit coupling, the one-body term can be written 

H {1) - V \h {1) ] d f d (2) 

L J no \ma 



n,7n,a 



where S na creates an electron of spin a on atomic orbital n of the molecule and [-H^oi] is a hermitian matrix. For 
simplicity, the atomic basis orbitals are taken to be orthonormal, so that the anticommutator {(V ncn d ma '} = 5 nm 8 aa i . 
Extension to non-orthogonal bases is straightforward in principle. For a general spin-rotation invariant two-body (e.g. 
Coulomb) interaction, 

H mr ^ — — ^ ^ U nm p n p m: (3) 



* mol 



n , rn 



where p n = Y, a d\ la d na . 

Due to their large density of states (and consequently good screening) , the macroscopic metallic electrodes (labeled 
a € [1, . . . , M}) may be modeled as non- interacting Fermi gases: 



M 

H loads = ^2 ^2 e *ff C L C *<" ( 4 ) 
a 

where ct creates an electron of energy efe CT in lead a. The electrostatic interaction of molecule and electrodes due 

to the electric dipoles formed at each molecule-electrode interface may be included in as discussed in Sec. IIIII 

Tunneling of electrons between the molecule and the electrodes is described by the Hamiltonian 

M 

^T = EEE (V nk di a c ka + H.c.) . (5) 

a=l tgn n,rr 



B. Non-equilibrium Green's functions 

The electronic system (jTJ) of molecule plus electrodes has an infinite Hilbert space. A formal simplification of the 
problem is obtained through the use of the Green's function a 23 i 24 

g„ ct ,w(<) = -tf(t)({<M*),dL'(°)}>. ( 6 ) 

G< a , m At) = iidlA^dnait)), (7) 

known as the retarded and Keldysh "lesser" Green's functions, respectively. Physical observables in the molecular 
transport junction can be expressed in terms of G and G < . In this section, an exact Dyson equation for G is derived, 
along with a corresponding Keldysh equation for G < . 

Setting h = 1, the retarded Green's function obeys the equation of motion^ 

i^-G n<J ,m*>{t) = 5{t)8 nm 5 aa , ~ iO{t) ({[d„ CT (t),ff junc tion] ,<*L'( )}) ■ ( 8 ) 
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For the purposes of this article, only the spin-diagonal term G na ^ ma is needed. Evaluating the commutators [d na , H mo i] 
and [d na , Ht], and noting that [d na , H\ ca ^ s ] = 0, Eq. ([5]) becomes 



d 

i- 



G„{t) = S(t)S nm + \h ( 2i] Gn'a, m «(t) + E E V <* 9k<r,m<r(t) + E ^'^i^W. ( 9 ) 

* — ■* l J n(j ,7i (T * — ' r — — 



where 

= -i6{t) ({c fcCT (t),4 CT ,(0)}) (10) 



and 



= "WW ({pn'(*)d™r(*),4ur/(0)}) • (H) 



The equation of motion for gka,ma(f) is 

i^- t 9k*, m *(t) = -»0(t) ({[ck*(t),H iunction ] • (12) 

Evaluating the commutators [c^ifr] and [cfe CT , -ffieadsL ancl noting that [c^o-, £f mo i] = 0, Eq. (I12p becomes 

i-Q^9ku,mcr(t) = i-ka 9ka ;mcr(t) + E ^n' fe^n'<7,m<7 (*) • (13) 

Fourier transforming Eqs. (J9j) and (fT3|) into the energy domain, and eliminating Qka^ma 

(E), one arrives at the 

following matrix equation for G(E): 

IS - #moi - St (E)] G(E) = 1 + UG^(E), (14) 
where the retarded tunneling self-energy matrix is 

FrWU-^-^EE d^F - (15) 

Eq. (fT^|) may be recast in the form of Dyson's equation 



G(E) = [IE - ffW _ St(£) _ Sc( ^)j ~ (16) 

via the ansatz UG^(E) = T*c(E)G(E), which defines— the retarded Coulomb self-energy matrix Y<c(E), the central 
quantity of the many-body theory, which must be determined via an appropriate approximation, as discussed below. 

For nonequilibrium problems, the Keldysh "lesser" self-energy and Green's function are also needed. G < is deter- 
mined by the Keldysh equation^! 

G<(E) = G(E)-£<(E)Gi(E), (17) 

where Y, < (E) — S^(_E) + T^^,(E), and G^{E) is the Fourier transform of the advanced Green's function. The "lesser" 
tunneling self-energy is 

[ S t(£)Uw = <£/«(^^UC*w. (18) 

a 

where f a (E) — {1 + exp[(i? — /ia)/ft.BT Q ]} _1 is the Fermi-Dirac distribution for lead a, and 

F a (E)] na , ma , = 2tt«W E KkVmk S(E - e ka ) (19) 

k£a 

is the tunneling- width matrix for lead a. 

Once the Green's functions are known, the relevant physical observables can be calculated. For example, the current 
flowing into lead a is given by2£ 

I a = j J dETr{T a (E)(G<(E) + f a (E)[G(E)-GHE)])}. (20) 
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FIG. 1: A schematic representation of detailed balance. The left-hand and right-hand sides of Eq. (|22[) are represented in (a) 
and (b), respectively. 



C. Sequential tunneling limit 

In the limit of infinitessimal lead-molecule coupling St/AibT — > 0, coherent superpositions of different energy 
eigenstates of the molecule can be neglected, and the junction Green's function becomes 

s hm o [G(E)] na<ma ee [G mol (£)W™ = E WW + W] E K tTX+iol (21) 

where \v) and \v') are many-body eigenstates of the isolated molecule satisfying -ff mo i|^} = E v \u), etc. The nonequi- 
librium probabilities V(v) can be determined by solving a system of semiclassical rate equations for sequential 
tunneling] 15 i 18 ' 19 i 32 For steady-state transport, they satisfy the principle of detailed balance (see Fig. [TJ: 

Viv) E f a"Va - ^) = W) E f t 1 " f° " ■ ( 22 ) 

Here the rate constants are given by Fermi's Golden rule as 

fW =T r{r Q (iv -E V )C U '}, (23) 

where T a (E) is given by Eq. (fl~9|) and 

[C^'W.W = (i/|d» ff |i/)(i/|4a'l I/ > C 2 ^) 

are many-body matrix elements i 33 ' 3 > 35 From the normalization of the many-body wavefunctions, the total resonance 
width Y vv = 5^„r„ I/ scales as ~ ^ Q Tr{r a }/iV, where AT is the number of atomic orbitals in the molecule. 
For strongly-correlated systems, there is an additional exponential suppression of Eq. (|24|) as N — > oo due to the 
orthogonality catastrophe. 35 ' 36 

Linear response transport is determined by the equilibrium Green's functions. In equilibrium, the solution of the 
set of Eqs. (f2"2")l reduces to 

V(u)=e-ft B "- ltN J/Z, (25) 

where Z is the grand partition function of the molecule at inverse temperature (3 and chemical potential fj,. 
Eq. (|2ip implicitly defines the Coulomb self-energy matrix £q in the sequential tunneling limit via 

(26) 

T,^ 1 is the high-temperature limit of the Coulomb self-energy (i.e., the limit /csT/Tr{r a } ^> 1), and describes 
intramolecular correlations and charge quantization effects (Coulomb blockade away from resonance). The nonper- 
turbative treatment of intramolecular correlations provided by exact diagonalization of H mo \ in Eq. (|2ip allows for an 
accurate description of the HOMO-LUMO gap, which is essential for a quantitative theory of transport. 



Gmol(E) — 



D. Correction to the Coulomb self-energy 

In general, the Coulomb self-energy matrix Sc = + AScj where ASc describes the change of the Coulomb 
self-energy due to lead-molecule coherence emerging at temperatures fc^T < Tr{r Q }. Using this decomposition of the 
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(a) (b) 



FIG. 2: The correction to the Coulomb self-energy AEc is given in the self-consistent Hartree-Fock approximation by the sum 
of (a) the Hartree (direct) term and (b) the Fock (exchange) term. The wavy line represents the Coulomb interaction, which 
has the form U nm S(r — r') on the Keldysh time contour. 



Coulomb self-energy, Dyson's equation (|I6p can be rewritten in the following useful form: 

G-\E) =G m 1 jE)- E T -AE C , (27) 

where the self-energy terms Et + AEc describe the effects of finite tunneling width. Et is given by Eq. (fTTI)) . Here 

we point out that AEc — unlike Sq — can be evaluated perturbatively using diagrammatic techniques on the Keldysh 
time-contour (cf. Fig. Such a perturbative approach is valid, in principle, at temperatures/bias voltages satisfying 
max{T, eV/ks} > Tk, where Tr- is the Kondo temperature^ 7 - — or when there is no unpaired electron on the molecule 
(such as within the HOMO-LUMO gap of conjugated organic molecules). 

Et + AEc m ay be thought of as the response of the junction to turning on the tunneling coupling T a . A sub- 
tlety in the perturbative evaluation of AEc is that the diagrams determining the Coulomb self-energy are typically 
formulated 3 -^ in terms of the Green's functions of the noninteracting system, G@)(E) = (IE — H^L — Et) -1 and 
= (IE — + £0 + ) _1 , while the response of the junction is determined by the full Green's function G(E) 
[cf. Eq. (|2"Tj)) ]. Evaluating the diagrams in Fig. fusing G^(E) and G^ ol (E) would yield a correction to the Coulomb 

self-energy with a pole structure unrelated to that of E^ , so that adding the two together would not yield a physically 
meaningful result. Our strategy is thus to calculate AEc by reformulating the terms in the perturbative expansion 
in terms of the full Green's function G(E) via appropriate resummations. This procedure is in general nontrivial, but 
the result for the Hartree-Fock correction is given in Sec. HID 21 below, based on physical arguments. Higher-order 
self-energy diagrams can be included in a similar fashion. Electron-phonon couplin g 25 i 38 ' 39 i 40 i 41 ' 42 can also be included 
in the many-electron theory via the self-energy terms E _ p h + AE C ph , where E e _ p h is given by the usual self-energy 
diagramS ) 25 i 39 ' 40 i 41 ' 42 and AE c ~ ph is the corresponding correction to the Coulomb self-energy. 



1. Elastic cotunnehng approximation: E = Et + E c 

Far from transmission resonances and for T 3> Tk, AEc can be neglected. This is the limit of elastic 
cotunnehng The Green's functions are given by Eqs. (fTtj)) and (fTT)) . with Ec = Eq . Note that [E^^ does 
not make a finite contribution to Eq. (fT7|) when Et is finite. The full NEGF current expression (|20|) then reduces to 
the multi-terminal Biittiker formula*^ 

^ = -rE/ dET Pa (E)[f p (E)-f a (E)}, (28) 

where the transmission function is given by^ 

Tp a (E) = TY{T a (E)G(E)Y [ \E)G J{ (E)} . (29) 

The elastic cotunnehng approximation is a conserving approximation — current is conserved [cf. Eq. (|28[)] and the 
spectral function obeys the usual sum-rule. 
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FIG. 3: Self-consistent Hartree-Fock correction to the Coulomb self-energy matrix of a diatomic molecule versus lead chemical 
potential, shown in the vicinity of the HOMO resonance. Here fio = (shomo +£lumo) /2 is the chemical potential of the 
isolated molecule, T = 300K, and Ih = T 2 = 0.2eV (see Sec.|In|. 



2. Self-consistent Hartree-Fock correction: Ec = Eq ' + AEq F 
In the Hartree-Fock (HF) approximation, the retarded Coulomb self-energy matrix is given by 



HF] 

C J nc 



and [S§ F ] < = 0. In general, 



and 



(30) 



(d na d mcr /) 



i 

27T 



dE 1 [G < (£')] m(T ' ! „ (J 



lira 



) = 'Yj'P(v)(y\d'! nrT d m<Tl \v). 



(31) 



(32) 



The Hartree-Fock correction is then AS§ F 



vHF 



VHF 

na, 



-►0- 

123 



The Feynman diagrams representing this correction 
and ()30|) - ()32p yields a conserving approximation in 



are shown in Fig. [2] A self-consistent solution of Eqs 
which the junction current is again given by Eq. (|28|) . 

The use of the interacting Green's functions G and G mo i in the evaluation of the Hartree self-energy is clearly 
justified on physical grounds, since this yields the classical electrostatic potential due to the actual nonequilibrium 
charge distribution on the molecule. The direct (Hartree) and exchange (Fock) contributions to AEc must be treated 
on an equal footing in order to cancel the unphysical self-interaction, justifying the use of the same interacting Green's 
functions in the evaluation of the exchange self-energy. The diagrams of Fig. [2] would be quite complex if expressed 
in terms of noninteracting Green's functions, because both G and G mo i involve Sq , which includes all possible 
combinations of Coulomb lines and intramolecular propagators. 

The self-consistent Hartree-Fock correction for a diatomic molecule is shown in Fig. [3] The parameters were chosen 
so that the resonance width ~ J2 a Tr{r Q }/7V = 0.2eV, where N is the number of atomic orbitals in the molecule. The 
elements of the matrix ASq F are largest near a transmission resonance, but vanish on resonance. This behavior can be 
understood by considering the molecular correlation functions shown in Fig. 3] Inclusion of the tunneling self-energy 
without the corresponding correction to the Coulomb self-energy leads to a charge imbalance on the molecule near 
resonance. This in turn leads to a Hartree correction which tends to counteract the charge imbalance. A corresponding 
behavior is found for the exchange correction and off-diagonal correlation function. A non-self-consistent calculation 
would yield a much larger correction AE§ F , indicating the important role of screening. It should be pointed out that 
a treatment of screening in linear response is not adequate near resonance. 

As shown in Fig. the transmission peaks and nodes are not shifted by the self-consistent HF correction, and the 
transmission phase is not changed significantly. Transport properties [cf. Eqs. (|28|) and (|29|) ] are therefore qualitatively 
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FIG. 4: Equilibrium correlation functions near the HOMO resonance for a diatomic molecule (same parameters as in Fig. 
[3]). Top: Total molecular charge (N) in three different approximations: E = Eq (isolated molecule in the grand canonical 
ensemble); E = Eq +Et (elastic cotunneling approximation); and E = Eq' +Et + AE§ F (self-consistent HF correction). The 
charging step is broadened in the elastic cotunneling approximation, compared to that of the isolated molecule, and acquires 
a 'knee' near resonance. In the self-consistent HF approximation, the charge imbalance AN = (N) — {N)\z T ->o is strongly 
screened near resonance, leading to a charging step as steep as that of the isolated molecule, but asymptotically approaches the 
elastic cotunneling result away from resonance. Bottom: The correlation function (d\ tT d,2a) in the same three approximations. 



better described in the clastic cotunneling approximation than are the correlation functions, and the approximation 
is quantitatively accurate in the cotunneling regime |/i — /i rcs | > T, max{T, eV/k^} > 2V- 

The tendency toward charge quantization near resonance is significantly increased by the self-consistent HF cor- 
rection, as shown in Fig. [5] The steepness of the self-consistent charging step is limited only by thermal broadening. 
This result is consistent with previous theoretical studie o 47 ' 48 i 49 of Coulomb blockade in metal islands and quantum 
dots, but inconsistent with the behavior of the Anderson model,— where singular spin fluctuations modify this generic 
behavior. 



III. MOLECULAR HETERO JUNCTION MODEL 



Hetcrojunctions formed from 7r-conjugated molecules are of particular interest, both because of their rele- 
vance for nanotechnology^ and their extensive experimental characterization^^ A semi-empirical 7r-electron 
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FIG. 5: Transmission probability and phase calculated with and without the self-consistent HF correction. Note that the 
transmission peaks and nodes are not shifted by AE§ F , but the width of the transmission resonance is reduced asymmetrically. 
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Hamiltonian 6 - ' 50 ' 51 can be used to model the electronic degrees of freedom most relevant for transport: 

Hmol = ^ ^ ^ri^ nC j^n(7 ^ 4 (^nm^<f lcr ^mcr ~t~ H.C.) -)- ^ ^ ~ QnQm: (33) 
n,a n,m,cr n : m 

where d n(J creates an electron of spin a in the 7r-orbital of the nth carbon atom, e n is the atomic orbital energy, 
and t nm is the hopping matrix element between orbitals n and m. In the 7r-electron theory, the effect of different 
side-groups is included through shifts of the orbital energies e n . The effect of substituents (e.g., thiol groups) used to 
bond the leads to the molecule can be include d 52 ' 53 in the tunneling matrix elements V n k [cf. Eq. 
The effective charge operator for orbital n is&i^ 5 - 

Qn = d ™ d ™ ~ C n«V a /e - 1, (34) 

a a 

where C na is the capacitive coupling between orbital n and lead a, e is the electron charge, and V a is the voltage 
on lead a. The lead-molecule capacitances are elements of a full capacitance matrix C, which also includes the 
intra-molecular capacitances, determined by the relation C nm = e 2 [U~ 1 ] nm - The values C na are determined by the 
zero-sum rules required for gauge invariance 5 ^ 

M 

^ Cnm + ^ C na = 0, (35) 

m a— 1 

and by the geometry of the junction (e.g. C na inversely proportional to lead-orbital distance). 
The effective interaction energies for 7r-conjugated systems can be writte n 50 ' 51 

U nm = S nm U + (1 - S nm ) — U ° , (36) 

eyfl + a(R nm /A) 2 

where Uq is the on-site Coulomb repulsion, a = (J7o/14.397eV) 2 , and R nm is the distance between orbitals n and m. 
The phenomenological dielectric constant e accounts for screening due to both the cr-electrons and any environmental 
considerations, such as non-evaporated solvent^ With an appropriate choice of the parameters t nm , Uq, and e, the 
complete spectrum of electronic excitations up to 8-10eV of the molecules benzene, biphenyl, and trans-stilbene in the 
gas phase can be reproduced with high accuracy 5 ^ by exact diagonalization of Eq. (|33|) . An accurate description of 
excited states is essential to model transport far from equilibrium. Larger conjugated organic molecules can also be 
modeled 5 ^ via Eqs. ([53 1) and ([M ]) . 

(T-orbitals can also be included in Eq. Q33p as additional energy bands, and the resulting multi-band extended 
Hubbard model can be treated using the same NEGF formalism sketched above. Tunneling through the c-orbitals 
may be important in small molecules ; 56 ' 57 especially in cases where quantum interference leads to a transmission node 
in the 7r-electron system. 

The biggest uncertainty in modelling single-molecule heterojunctions is the lead-molecule coupling. 58 For this rea- 
son, we take the two most uncertain quantitites characterizing lead-molecule coupling — the tunneling width T and the 
chemical potential offset A/i of isolated molecule and metal electrodes — as phenomenological parameters to be deter- 
mined by fitting to experiment. In the broad-band limits for the metallic electrodes, and assuming each electrode is 
covalently bonded to a single carbon atom of the molecule, the tunneling- width matrix reduces to a single constant: 
r" m (£) = r a 5 na S ma , where a is the orbital connected to electrode a. Typical estimate o 52 ' 53 ' 60 indicate T < leV for 
organic molecules coupled to gold contacts via thiol groups. 

IV. BENZENE(l,4)DITHIOL JUNCTION 

As a first application of our many-body theory of molecular junction transport, we consider the benchmark system 
of benzene(l,4)dithiol (BDT) with two gold leads i 27 ' 28 ' 29 ' 61 ' 62 ' 63 ' 64 ' 65 The Hamiltonian parameters for benzene are^i 
Uq = 8.9eV, e = 1.28, t nm = 2.68eV for n and m nearest neighbors, and t nm = otherwise. We consider a symmetric 
junction (symmetric capacitive couplings and ri = T2 = T) at room temperature (T=300K). 

A. Linear electric and thermoelectric junction response 

Thermoelectric effect o 27 ' 28 ' 66 provide important insight into the transport mechanism in single-molecule junctions, 
but are particularly sensitive to correlations ; 67 ' 68 calling into question the applicability of single-particle or mean-field 
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FIG. 6: Thermoelectric power of a BDT-Au junction at T — 300K as a function of the lead-molecule chemical potential 
mismatch for three different tunneling widths. Comparison to the experimental valued 7.0±0.2/xV/K fixes /xa u = /xo — (3.22 ± 
0.04)eV (see inset showing closeup of experimentally relevant parameter range). 



theory. We can now investigate thermoelectric effects in molecular heterojunctions for the first time using many-body 
theory. The thcrmopower S of a molecular junction is obtained by measuring the voltage AV = —SAT created 
across an open junction in response to a temperature differential AT. In general, S can be calculated by taking the 
appropriate linear-response limit of Eq. (|20p. which includes both elastic and inelastic processes. However, for purely 
elastic transport [cf. Sees. Ill D II and III D 2] . Eq. j28]) can be used to derive the well-known result 

S(», T) = --==— V ' r . (37) 

eT sr^E) (-% )dE 



Thermopower measurements^^ provide a means to determine^ the lead-molecule chemical potential mismatch 
A^i = /! Au - A*o, where /x = (ehomo + £lumo) /2 and £homo (elumo) is the HOMO (LUMO) energy level. Figure 
|6] shows the thermopower of a BDT junction as a function of A/x, calculated from Eqs. ([29]) and (|37|) in the elastic 
cotunneling approximation. As pointed out by Paulsson and Datta^ S is nearly independent of T away from the 
transmission resonances, allowing A/x to be determined directly by comparison to experiment. The thermopower of a 
BDT-Au junction was recently measured by Baheti et a/.,— who obtained the result S = (7.0 ± 0.2)/xV/K. Equating 
this experimental value with the calculated thermopower shown in Fig.[6l we find that — 3.25eV < fJ.Au — fJ-0 < — 3.15eV 
over a broad range of T values. The Fermi level of gold thus lies about 1.8eV above the HOMO resonance, validating 
the notion that transport in these junctions is hole-dominated. Nonetheless, /xa u — £homo is sufficiently large that 
the elastic cotunneling approximation is well justified. 

The only other free parameter in the molecular junction model is the tunneling width T, which can be found by 
matching the linear- response conductance to experiment. Although there is a large range of experimental values j« the 
most reproducible and lowest resistance contacts were obtained by Xiao, Xu and Tao,— who report a single-molecule 
conductance value of 0.011G . This fixes T = (0.63 ± 0.02)eV with ^x Au - /x = -3.22 ± 0.04eV. This value of T 
is within the range predicted by other groups for similar molecules j 52 ' 53 i 58 ' 60 With the final parameter in the model 
fixed, we can now use our many-body transport theory to predict the linear and nonlinear response of this molecular 
junction. 

The transmission probability T{E)\e=i_i of the BDT junction is shown as a function of fj, — /xau in the upper 
panel of Fig. [7] This is the linear-response conductance in units of the conductance quantum Go = 2e 2 /h. The 
chemical potential /x is related to the gate voltage V g in a three-terminal junction via the gate capacitance C g . The 
transmission spectrum shown in Fig. [7] exhibits several striking features: a large but irregular peak spacing with 
an increased HOMO-LUMO gap, non-symmetric, Fano-like resonance lineshapes, and transmission nodes due to 
destructive interference in the coherent quantum transport. 
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u^ Au (eV) 

FIG. 7: Linear-response conductance versus chemical potential (top panel) and differential conductance versus chemical 
potential and bias voltage (bottom panel) for a benzene(l,4)dithiol-Au junction. The calculation was carried out within the 
elastic cotunneling approximation using Eq. (|28[) , with I\ = F2 = 0.6eV and T = 300K. Here Go = 2e 2 /h is the conductance 
quantum. Note the asymmetric Fano-like lineshapes and the transmission nodes (arising from destructive interference) in the 
top panel. The numbers within the central diamonds of the bottom panel indicate the quantized charge on the molecule 
(relative to neutral benzene). The full spectrum is shown for completeness. 



B. Non-linear junction transport 

We next calculate the differential conductance dl/dVbias of the BDT-Au junction as a function of [i and Vbias 
(see Fig. [3 lower panel). The current was calculated in the elastic cotunneling approximation using Eqs. I]28p and 
(1291) . This approximation accurately describes nonresonant transport, including the transmission nodes, as well as the 
positions and heights of the transmission resonances, as discussed in Sec. lIID~2l The differential conductance spectrum 
of the junction exhibits clear signatures of excited-state transpor t 30 i 71 and an irregular "molecular diamond" structure 
analogous to the regular Coulomb diamonds observed in quantum dot transport experiments^ The charge on the 
molecule is quantized within the central diamonds of Fig. an important interaction effect inaccessible to so-called 
ab initio mean-field calculations. In Fig. the full spectrum is shown for completeness, although the junction may 
not be stable over the entire range of bias and gate voltages shown. 

Apart from the central HOMO-LUMO gap, the widths of the diamonds in Fig. [7] can be roughly explained via a 
capacitive model in which the molecule is characterized by a single capacitance C mo i = e 2 / (U nm ), where (U nm ) — 
5.11eV is the average over all molecular sites. The HOMO-LUMO gap of ^10eV is significantly larger than this 
estimate of the charging energy, an indication of the significant deviations of our theory from a simple constant 
interaction model. 

Charge quantization, also known as Coulomb blockade, has been observed in several different types of molecular 
heterojunctionS ) 72 ' 73 i 74 ' 75 i 76 but has not yet been observed in BDT junctions due to the difficulty of gating such 
small molecules^ The unambiguous observation of Coulomb blockade in junctions involving larger molecules (with 
smaller charging energies) indicates that such interaction effects, lying outside the scope of mean-field approaches, 
are undoubtedly even more pronounced in small molecules like BDT. Important aspects of this phenomenon remain 
to be understood in larger molecules, such as the anomalously low reported values of the charging energyi 72 i 73 i 74 i 75 i 76 

A zero-bias cross section of the bottom panel of Fig. [7] reproduces the transmission spectrum shown in the top 
panel of the same figure. Increasing the bias voltage, we find that the resonances split into negatively sloped particle- 
like (\e\dV/dn = - (Ci + C 2 + C g ) /d) lines and positively sloped hole-like (\e\dV/dfi = + {C\ + C 2 + C g ) /C 2 ) 
lines, where G\ and C 2 are the mean lead-molecule capacitances, defined by C a — (C na ), a = 1,2. In the symmetric 
coupling case (Ci = C 2 with C g jC\ <C 1) the lines therefore have slopes of -2 and +2 for particle-like and hole-like lines, 
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respectively. Within the V-shaped outline traced by the particle-like and hole-like lines, we find signatures of resonant 
tunneling through electronic excited states in the many narrow, nearly parallel resonance lines. While transport 
through electronic excited states has not yet been unambiguously identified in single-molecule heterojunctions, it has 
been observed in quantum dots 30 -^ and carbon nanotubes^ 

An accurate description of the HOMO-LUMO gap is essential for a quantitative theory of transport in molecular 
heterojunctions. The central HOMO-LUMO gap shown in Fig.[7]is significantly larger than that predicted 1 - by density 
functional theory — which neglects charge-quantization effects — but is consistent with previous many-body calculations 
in the sequential tunneling regime i 15 ' 18 ? 19 It should be emphasized that the transport gap A/zt r in a molecular junction 
exceeds the optical gap ?iw m i n of an isolated molecule. Roughly speaking, A/i tr ~ huj m i n + e 2 /C mo \ + E x , where e 2 /C mo i 
is the charging energy (see discussion above) and — E x is the exciton binding energy. Excitonic states of the BDT-Au 
junction can be identified in the differential conductance spectrum of Fig. as the lowest-energy (i.e., smallest bias) 
excitations outside the central diamond of the HOMO-LUMO gap, from which it is apparent that A/i tr ~ 2?kj m i n for 
a BDT-Au junction. 

V. CONCLUSIONS 

In conclusion, we have developed a many-body theory of electron transport in single-molecule heterojunctions that 
treats coherent quantum effects and Coulomb interactions on an equal footing. As a first application of our theory, 
we have investigated the thermoelectric power and differential conductance of a prototypical single-molecule junction, 
benzencdithiol with gold electrodes. 

Our results reproduce the key features of both the coherent and Coulomb blockade transport regimes: Quantum 
interference effects, such as the transmission nodes predicted within mean-field theory ^ 57 i 79 are confirmed, while 
the differential conductance spectrum exhibits characteristic charge quantization "diamonds" 72 i 73 i 74 i 75 i 76 — an effect 
outside the scope of mean-field approaches based on density-functional theory. The HOMO-LUMO transport gap 
obtained is consistent with previous many-body treatments in the sequential tunneling limit i 15 i 18 ' 19 

The central object of the many-body theory is the Coulomb self-energy Sc of the junction, which may be expressed 
as Sc = + ASci where is the result in the sequential tunneling limit, and ASq is the correction due to a 

finite tunneling width F. In this article, we have evaluated E^? exactly, thereby including intramolecular correlations 
at a nonperturbative level, while the direct and exchange contributions to AEc were evaluated self-consistently using 
a conserving approximation based on diagrammatic perturbation theory on the Keldysh contour. An important 
feature of our theory is that this approximation for AEc can be systematically improved by including additional 
processes diagrammatically. In this way, important effects such as dynamical screening, spin-flip scattering^ 7 , and 
electron-phonon couplin g 25 ' 38 i 39 ' 40 i 41 ' 42 can be included as natural extensions of the theory. 

Our general theory of molecular junction transport should be contrasted to results for transport in the Anderson 
mode h ' The Anderson model provides important insight into nonperturbative interaction effects in quantum 
transport through nanostructures; however, the internal structure of the molecule is neglected. Moreover, since it 
is limited to a single spin-degenerate level, it can only describe a single Coulomb diamond with an odd number of 
electrons on the molecule, and is therefore not applicable to transport within the HOMO-LUMO gap of conjugated 
organic molecules, where the number of electrons is even (more precisely, the HOMO-LUMO gap is taken to be infinite 
in the Anderson model). 

Recently, a NEGF approach for molecular junction transport was introduced 8 -^ using Green's functions that are 
matrices in a basis of many-body molecular eigenstate outer products it is unclear what relation this approach 

has to the approach proposed herein involving single-particle Green's functions, which are matrices in the atomic 
orbital basis. 
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